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ABSTRACT 

We estimate the luminosity evolution and formation rate for over 900 GRBs 
by using redshift and luminosity data calculated by Band, Norris, & Bonnell 
(2004) via the lag-luminosity correlation. By applying maximum likelihood tech- 
niques, we are able to infer the true distribution of the parent GRB population's 
luminosity function and density distributions in a way that accounts for detec- 
tor selection effects. We find that after accounting for data truncation, there 
still exists a significant correlation between the average luminosity and redshift, 
indicating that distant GRBs are on average more luminous than nearby coun- 
terparts. This is consistent with previous studies showing strong source evolution 
and also recent observations of under luminous nearby GRBs. We find no evi- 
dence for beaming angle evolution in the current sample of GRBs with known 
redshift, suggesting that this increase in luminosity can not be due to an evolu- 
tion of the collimation of gamma-ray emission. The resulting luminosity function 
is well fit with a single power law of index Z/ -1,5 , which is intermediate between 
the values predicted by the power-law and Gaussian structured jet models. We 
also find that the GRB comoving rate density rises steeply with a broad peak 
between 1 < z < 2 followed by a steady decline above z > 3. This rate density 
qualitatively matches the current estimates of the cosmic star formation rate, fa- 
voring a short lived massive star progenitor model, or a binary model with a short 
delay between the formation of the compact object and the eventual merger. 

Subject headings: gamma rays: bursts — gamma rays: theory 
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1. Introduction 

There are currently roughly three dozen gamma-ray burst events (GRBs) for which we 
have independently measured redshifts. Most of these redshift determinations come from 
either identification of absorption lines in the afterglow spectra, attributed to the gas in the 
host galaxy, or from observations of emission lines from the host galaxy. The combination 
of these techniques has resulted in a small but growing GRB sample with redshifts ranging 
from 0.0085 to 4.5 and a distribution peaking around z ~ 1. From this small sample, it is 
already abundantly clear that the isotropic equivalent energy E iso released in the prompt 
GRB phase is not a standard candle. The total radiated energy taken at face value (i.e. 
when not correcting for a beaming factor dfl) clearly spans several orders of magnitude, 
ranging from 10 47 for the closest event, GRB 980425 at z = 0.0085 (Kulkarni et al. 1998), to 
10 54 for GRB 990123 at z = 1.6004 (Kulkarni et al. 1999). Recently Sazonov, Lutovinov, & 
Sunyaev (2004) and Soderberg (2004) have reported on gamma-ray observations of a nearby 
underluminous GRB occurring at redshift of z = 0.106. These new findings have added to 
the speculation that there is either a substantial under luminous population of GRBs which 
cannot be seen at large distances and/or that nearby events (z < 0.15) are underluminous 
compared to distant counterparts, pointing to the evolution of the average energy emitted 
by a GRB with time. 

A measure to the extend to which luminosity evolution exists in the GRB population, 
along with their true luminosity function and density distribution, may yield important clues 
regarding the nature of gamma-ray bursts and how they're progenitors have evolved with 
time. Although the physics of the underlying GRB engine is hidden from direct observa- 
tion and is yet uncertain, the total GRB energy budget is most likely linked to the mass 
and/or rotational energy of the GRB progenitor. Understanding of how this energy budget 
has changed with time may offer constraints on progenitor properties and may ultimately 
point to the physics leading to their explosions. Since GRB progenitors are most likely 
linked to compact objects (supermassive rotating star, black hole or neutron star mergers) 
understanding how the GRB luminosity function evolves with time may give insight to the 
host environment in the early universe, namely the star formation rate or initial stellar mass 
functions at high redshifts. 

Any attempt at quantifying the evolution of intrinsic source properties of parent pop- 
ulations must account for Malquist type biases. Detection thresholds prevent events below 
a certain flux from being observed, resulting in the detection of only bright objects at large 
distances. Combined with the fact that bright events are typically rare, it is very easy in as- 
tronomy to incorrectly conclude that the distant universe is filled with extremely bright rare 
objects. Any attempt at measuring the correlation between luminosity and redshift without 
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properly accounting for selection effects will grossly overestimate the correlation strength 
between the two variables. Flux limited samples are a classic problem in astronomy, which 
manifested prominently in early quasar studies. Fortunately, straight forward methods have 
been devised to account for such effects based on maximum likelihood techniques. These 
methods allow for the correct estimation of the correlation strength between a truncated 
data set as well as an estimate on the underlying parent population. The "catch" of such 
techniques is that the overall normalization of the resulting parent distributions cannot be 
determined, although their functional forms are constructed in such a way to account for the 
data truncations. 

These techniques also have to limitation of requiring a large sample sizes and more 
importantly, an extremely good understanding of the survey's detection thresholds (i.e. the 
flux cutoff for magnitude limited samples). The use of the current sample of GRBs with 
known redshift is limited by both of these restrictions. The current size of a little over two 
dozen bursts does not lend itself well to producing statistically robust results, especially in 
the high and low redshift regimes for which only a handful of events have been detected. 
Furthermore, the sample is an accumulation of observations from several different spacecraft, 
all of varying detector thresholds. It would seem that these limitations could only be over- 
come by the accumulation of a larger data set with consistent detector thresholds which is 
expected to come from the Swift spacecraft and the upcoming GLAST mission. 

Fortunately, several authors have announced empirical Cepheid like correlations linking 
intrinsic burst properties, such as luminosity (Norris, Marani, & Bonnell 2000) and the total 
radiated energy (Amati et al. 2002; Ghirlanda et al. 2004a) to other GRB observable. These 
correlations may allow for the determination of burst redshifts directly from the gamma-ray 
data, which has the advantage of being relatively insensitive to extinction and observable at 
far greater distances than afterglow line measurements. The first of these correlations was 
reported by Norris, Marani, & Bonnell (2000). Using 6 BATSE detected bursts with known 
redshift, they found an anti-correlation between the source frame lag between the 25-50 keV 
and 100-300 keV emission and the absolute luminosity of the GRB. More recently, (Ghirlanda 
et al. 2004a) reported an empirical correlation between the collimation correction total energy 
E 1 radiated by the burst and the rest frame energy at which most of the prompt radiation is 
emitted E pk . Using these relationships, it is now possible to estimate "pseudo" redshifts for 
a much larger number of GRBs detected by the BATSE instrument which perviously lacked 
any information as to their distance. More importantly, the BATSE detector threshold is 
relatively well understood for the entire sample, making the resulting pseudo redshift data 
excellent for statistical analysis. 

In this paper, we examine the issue of luminosity and density evolution by using a 
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sample of over 900 BATSE GRBs for which the luminosity and redshift where recently esti- 
mated by Band, Norris, Bonnell (2004) through the use of the lag-luminosity correlation. We 
limit our analysis to the lag-luminosity correlation primarily due to the lack of jet opening 
angle 0j information that is required for the use of the E^-E p k relation. This relationship 
requires knowledge of 0j in order to determine the collimation factor, which is only known 
for bursts with measured jet break times and hence cannot be used with the BATSE sam- 
ple in consideration for this paper. We found that the more general, and much broader, 
correlation between intrinsic E pk and E iso reported by Amati et al. (2002) did not provide 
redshift constraints for a majority of the bursts in our sample. This is consistent with recent 
observations by Nakar & Piran (2005) and Band & Preece (2005) who also found large 
fractions of their BATSE samples to be inconsistent with the Amati correlation. Therefore 
we limit the current analysis to distances estimated through the use of the lag-luminosity 
relationship. 

To our sample, we apply statistical techniques developed by Lynden-Bell (1971) and 
Efron & Petrosian (1992) and first applied to GRB analysis by Lloyd-Ronning, Fryer, & 
Ramirez-Ruiz (2002) to measure the underlying luminosity and density distribution in a 
way that properly accounts for the detection thresholds of the BATSE instrument. We find 
a strong (11.63 a) correlation between luminosity and redshift that can be parameterized 
as L(z) = (1 + z) 1 - 7 ^- 3 . The resulting cumulative luminosity function N(L') is well fit by 
double power law separated by a break energy of about 10 52 ergs s™ 1 , with the differential 
luminosity function dN/dL' exhibiting a power law shape of L~ L5 below this luminosity. 
We show that the GRB comoving rate density increases roughly as Pj(z) oc (1 + z) 2 - 5 to a 
redshift of z ~ 1 followed by a flattening and eventual decline above z > 3. This rate density 
is in qualitatively agreement with recent photometric estimates of the cosmic star formation 
rate (SFR), as would be expected from massive short lived progenitors. 

In §2, we describe the data set that we use in our study. In §3 we discuss the statistical 
methods applied to this data to estimate the GRB luminosity function and comoving rate 
density as well as to test for any correlation between luminosity and redshift. In §4 we 
present the resulting demographic distribution functions of this analysis followed in §5 by 
a discussion of the implications of the shape and evolution of the luminosity function and 
comoving rate density on various jet profile. We show that there is no evidence for beaming 
angle evolution in the current sample of GRBs with known redshift, suggesting that the 
variation of the observed luminosity with redshift can not be due to an evolution of the 
collimation of gamma-ray emission. We conclude by examining how the similarity between 
the SFR and the GRB comoving rate density tentatively favors short lived progenitor models. 



- 5 - 



2. Data 

For this analysis we utilize data for 1438 BATSE detected GRBs presented in Band, 
Norris, Bonnell (2004), hereafter BNB04. This sample includes peak photon flux f p k in the 
50-300 keV band on a 256 ms timescale, the burst duration T 90 , and measured lags and 
their uncertainties for each burst. From these lag measurements, the authors infer each 
burst's luminosity and redshift by use of the lag-luminosity correlation, allowing also for an 
estimation of the intrinsic E pk and E iso for each burst. Of these 1438 bursts, 1218 have 
positive lags making them suitable for this analysis. This data is shown in Figure 6 with an 
imposed flux cut set at 0.5 photon cm -2 s -1 , leaving a total of 985 bursts. 

The lags measurements used in this sample where made using a cross-correlation analy- 
sis similar to that previously employed by Band et al. (1993) and Norris, Marani, & Bonnell 
(2000). The cross correlation method has been widely used in x-ray and gamma-ray as- 
tronomy, and is well suited for timing analysis between two signals. In this application, the 
normalized discrete cross correlation function is given by 

CCFjr) = yM) (i) 

where t' is commonly referred to as the lag between fit) and git) and aj = {fit) 2 ) 1 / 2 . 
By maximizing the CCF function (i.e. by maximizing the area of the product of the two 
functions) as a function of t', an estimate of the temporal offset of the two signals can be 
made. If g(t) leads f(t) by to (i-e. f(t) = git + to)) than the CCF curve peaks at t' = t . 

In BNB04, the authors utilize 64ms count data gathered by BATSE's Large Area Detec- 
tors (LADs) which provide discriminator rates with 64 ms resolution from 2.048 s before the 
burst to several minutes after the trigger (Fishman et al. 1994). The discriminator rates are 
gathered in four broad energy channels covering approximately 25-50, 50-100, 100-300, and 
300 to about 1800 keV allowing for excellent count statistics since the photons are collected 
over a wide energy band. BNB04 measure the temporal offset or lag between channel 3 
(100-300 keV) and channel 1 (25-50 keV) light curves to produce the CCF31 lags listen in 
their sample. 

The shifting of the GRB spectra out of (or into) the observers frame, otherwise known 
as the k-correction, was accounted for in the analysis performed in BNB04. They perform 
spectral fits for most of the bursts in their sample and for those which cannot yeild a fit, a 
"Band" spectral model with average parameters is assumed for the spectra. The effects of 
time dilation and k-correction are then used to obtain the source frame lag and also applied 
to the energy flux to obtain a bolometric luminosity. 
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As was the case in the original Norris, Marani, & Bonnell (2000) paper, the CCF method 
used in BNB04 can result in lag measurements which are less than the 64ms time resolution 
of the BATSE instrument. In these cases, the associated errors of these values tend to 
be quite large, reducing the significance of their associated luminosity and redshift values. 
These errors are taken into consideration in the maximum likelihood techniques performed in 
our analysis. Therefore, bursts with extremely short lags (and hence high luminosity's) are 
weighted accordingly. A plot of the lag-luminosity plane for the events under consideration 
along with the errors in the lag measurements are shown in Figure 6. 



2.1. Estimating Redshifts 

Using these lag measurements, BNB04 utilize the lag-luminosity correlation to estimate 
the luminosity of each event. This empirical correlations was reported by Norris, Marani, & 
Bonnell (2000) who used the CCF method to measure the lag between BATSE's channel 3 
and channel 1 energy light curves for 6 GRBs with independently measured redshift. They 
concluded that there was an anti-correlation between the source delay in the low and high 
energy emission and the absolute luminosity of the GRB showing that high luminosity events 
exhibited very small intrinsic (source frame) lag, whereas fainter bursts exhibited the largest 
time delay. This empirical correlation can be expressed as 

L = 2.51 x 10 51 (A*70.1) -1 - 15 (2) 

where At' is the source frame lag related to the observed lag At' obs by a time dilation factor 
of (1 + z)~ l . The fact that the lag-luminosity correlation relates two source frame quantities 
(i.e. luminosity and intrinsic lag) would make it seem that knowledge of the redshift is needed 
a priori. As it turns out this is not the case. A simple numerical iteration routine can be 
used to solve for the redshift of a GRB which lacks any information as to its distance. This 
is done by first making an initial guess for z (say z ~ 1) to obtain the lag in the comoving 
frame At' = At' o6s /(l + z). This in turn gives us an initial value for the luminosity through 
the use of the lag-luminosity relation. This luminosity is then used in combination with the 
burst's energy flux to obtain a value for the luminosity distance Dl through the standard 
relation 

D L - M (3) 

V J256 

where /256 is the peak flux in the 256 ms timescale and dQ is the beaming factor. This 
distance is then compared to the Dl that can be calculated directly from the guessed redshift 
z by assuming standard cosmological parameters (H = 65 km s _1 , Vt m = 0.3, Q\ = 0.7) 
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and using the expression 



= + / 

-no Jo 



dz 



(4) 



The value for z is then varied until the luminosity distances obtained from the two separate 
methods converge to within some predetermined precision. 

We note that it has been suggested by Salmonson (2001) and Norris (2002) that the 
lag-luminosity relationship should be a broken power law in order to accommodate GRB 
980425. This burst was associated with SN 1998bw and when using the distance to the 
supernova, the GRB appears under luminous compared to the other bursts that fall on the 
lag-luminosity correlation. In their analysis, BNB04 note that this break has been suggested 
to fit a single point, which may or may not be associated with the SNe event and hence 
decide to use a single power law of -1.15. 

The physical origin of the lag-luminosity correlation is not immediately clear. Funda- 
mentally, this observed lag is due to the evolution of the GRB spectra to lower energies, so a 
relationship between the rate of spectral decay and luminosity is expected Kocevski & Liang 
(2003). This implies that the mechanisms resulting in the "cooling" of the GRB spectra is 
intimately related to the total energy budget of a GRB or its collimation factor. Other pur- 
posed theories attempt to explain the lag- luminosity correlation as being due to the effect of 
the viewing angle of the GRB jet (Kobayashi, Ryde, & MacFadyen 2002; Ioka & Nakamura 
2001), and or kinematic effects (Salmonson 2000). In any case, the use of this correlation is 
similar to methods used to calibrate Type la supernova luminosities based on the empirical 
correlation between their peak magnitude and rate of light curve decay (Phillips 1999). The 
lack of a clear physical interpretation of these correlations does not immediately preclude 
their use in determining, or refining, luminosity estimates. 



The luminosity and redshift data calculated by BNB04 gives us an enormous sample from 
which to investigate the evolution of the GRB luminosity function. As with any cosmological 
source, it is important and revealing to understand of how the average luminosity and density 
has evolved with cosmic time. Attempting to do so by simply measuring the correlation 
coefficient between the flux truncated luminosity and redshift data in the BNB04 sample 
without properly accounting for the detector selection effects would grossly overestimate 
the correlation strength. This is true whenever an estimate of correlation is made between 
two variables that suffer from data truncations, with the resulting correlation coefficient 
representing the truncation itself and not the underlying relation. 



3. Analysis 
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There have been several methods developed in astronomy to account for such selection 
effects, based largely on maximum likelyhood techniques (see Petrosian 1992 for a review). 
In our analysis, we use a nonparametric statistical technique originally proposed by Lynden- 
Bell (1971) for applications in flux limited quasar studies. This so called C-Method has 
been used successfully to reconstruct underlying parent distributions for quasars and GRBs 
samples by Maloney & Petrosian (1999), and Lloyd- Ronning, Fryer, & Ramirez- Ruiz (2002) 
respectively. The parent luminosity and redshift distributions which the method estimates 
allows for the construction of a GRB luminosity function, a measure of the number of bursts 
per unit luminosity, and an estimate on the comoving rate density, a measure of the number 
of bursts per unit comoving volume and time. 

The C-Method has two important limitations, or stipulations, to its use. First, the 
truncation limit below which no observations can be made must be well known. This is 
not a problem in our case, since the detector threshold of the BATSE instrument is well 
understood and BNB04 quantify the truncation limit of their sample. Secondly, the parent 
luminosity and redshift distributions can only be estimated in a bivariate manner if the two 
variables are uncorrelated. This is a limitation of all nonparametric techniques which rely on 
the assumption of stochastic independence. Therefore, it is necessary to first determine the 
degree of correlation between the two variables, in our case luminosity and Z — 1 + z, and 
then produce an uncorrelated data set through the transformation L — > L' = L/g(z), where 
g(z) parameterizes the luminosity evolution. Using this uncorrelated data set, it is then 
possible to apply the C-Method to estimate the underlying parent luminosity and redshift 
distributions. To estimate the degree of correlation we use a simple test of independence for 
truncated data put forth by Efron & Petrosian (1992) which is based in part on Lynden- 
Bell's C-Method. Below we describe the details of both Lynden-Bell's C-Method and the 
Efron & Petrosian independence test and how they are applied in our analysis. 



3.1. Test of Independence 

If the variables x and y in a data set are independent, then the rank Ri of X{ within 
that set should be distributed uniformly between 1 and iV with an expected mean E = 
(1/2)(JV + 1) and variance V = (l/T2)(iV 2 — 1). It is common practice to normalize the 
rank Ri such that for independent variables Ri has a mean of and a variance of 1 by 
defining the statistic Tj = (Ri — E) /V . A specialized version of the Kendell r statistic can 
be constructed to produce a single parameter whose value directly rejects or accepts the 
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hypothesis of independence. This quantity is commonly denned as 

_ — E) 



(5) 



Based on this definition, a r of 1 indicates ale correlation whereas a r of signifies 
a completely random data set. See Efron & Petrosian (1992) for a more detailed (and 
elucidating) proof of the applicability of normalized rank statistics. 

The modified version of this rank statistic proposed by Efron & Petrosian (1992) to 
test the independence of truncated data is based on a simple concept. Instead of measuring 
the ranks Ri for the entire set of observables, rather deal with data subsets which can be 
constructed to be independent of the truncation limit suffered by the entire sample. This is 
done by creating "associated sets" which include all objects that could have been observed 
given a certain limiting luminosity. We can define an associated set as 

Ji = {j : Lj > Li, Lu m j < Li} (6) 

In other words, for each burst % there can be constructed a data subset that includes all 
events within the range Li < L < oo and < z < z max (Li). The boundaries of an associated 
set for a given burst % are shown as dotted lines in Figure 6. In this scenario, we expect the 
rank Ri of Zi within the associated set 

Ri = {j E Ji : Zj < z^ (7) 

to be uniformly distributed between 1 and Nj, where Nj is the number of points in the 
associated set J«. Using these new ranks, we can again construct the mean and variance, 
except that now we replace N with Nj such that E = (1/2) (JV,- + 1) and V = (l/12)(iVj - 1). 
The specialized version of Kendell's r statistic is now given by 

r = (8) 

where the mean and variance are calculated separately for each associated set and summed 
accordingly to produce a single value for r. This parameter represents the degree of correla- 
tion for the entire sample with proper accounting for the data truncation. With this statistic 
in place, it is a simple matter to find the parametrization that best describes the luminosity 
evolution. This is accomplished by first choosing a functional form for the luminosity evolu- 
tion, which in our case we choose a simple power law dependence g(z) = (1 + z) a . We can 
then make the transformation L — > U — L/g{z) and vary a until r — > 0. 

An example of how well these methods are able to estimate underlying correlations in 
truncated data is shown in Figure 6. Here we have plotted a distribution of fake luminosity 
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and redshift data with some known power law dependence L oc (1 + z) p which is subjected 
to a flux cut L Hm oc (1 + z) q represented by the red dashed line. The crosses show the 
observable data whereas the dots represent the data that would otherwise be undetectable. 
The long dashed line is the best fit to the truncated data without any knowledge of the 
flux cut whereas the dash dot line is the reconstructed correlation when taking into account 
the flux threshold. This method fails when the undetected data points become significantly 
larger than the observable data set, with the exact boundary at which this occurs depending 
on the difference in the power law indices between the underlying correlation and the flux 
threshold. Since these quantities cannot be known a priori, it is explicitly assumed that a 
large data sample contains a sufficient amount of events above the flux threshold for the 
method to work. A histogram of the difference between the known correlation index and the 
reconstructed index (p — q) for multiple such simulations is shown in Figure 6. The error, or 
difference between the known p and the measured q is peaked about zero with a fwhm which 
roughly matches that error estimates that correspond to the 1 a range for this parameter 
given by the condition |r| < 1. 



3.2. Determination of Distribution Functions 

Once a parametric form that removes the the correlation between L and z is known, it 
is possible to use nonparametric maximum likelyhood techniques to estimate the underlying 
parent luminosity and redshift distributions. This luminosity distribution represents 
the cumulative GRB luminosity function with the redshift distribution cr(zi) representing the 
GRB density evolution. Petrosian (1992) has shown that many, if not most, of the familiar 
nonparametric methods used in astronomy to produce and cr(zi) reduce fundamentally 

to Lynden-Bell's C-Method. Consider the area, or number of events, in the box produced 
by the associated set shown in Figure 6. If Aq represents the number of points with L > L±, 
then let dN\ represent the number of points in the infinitesimal column between L\ and 
Li + dL\. The general premise behind the C-Method is that if the two variables (L, z) 
are stochastically independent, then the ratio between Aq and dAq should equal the ratio 
between <i<3> and the true cumulative distribution function <E>(Li) 

dNi d$ , . 

*r=*r < 9 > 

which can then be integrated to find <&{L). In the case of discrete data points, this integration 
becomes a summation, yielding the solution 

*m = n ( x + 4) ( iq ) 
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where Nj, is the number of bursts in the box defined by < z < z max (Lj) and Lj < L < oo. 
The value Nj is the same as Lynden-Bell's Cj in that it does not count the Lj object that 
is used to form the associated set. Similarly, we can construct the underlying cumulative 
redshift distribution function a{zj) by reversing the definition of the associated set such that 
Mj represents the number of bursts in the box < z < Zi and L min (zi) < L < oo. Then 

= ri 

k=2 

As mentioned in § 1, there are several important limitations to the C- method. First, 
the overall normalization of and cr(zi) is arbitrary, so information regarding the ab- 

solute numbers and densities cannot be obtained. Despite this, the shape of the bivariate 
distribution is constructed in such a way that it accounts for the data truncations. Due to 
this limitation, all distributions presented in this paper will have arbitrary normalizations. 
Secondly, it is clear from Equation 10 and 11 that the cumulative distribution function is 
not defined when either Nj or Mj are zero. This limitation restricts the use of the C-method 
to samples with a data size sufficiently large to ensure that all associated sets greater than 
j — 2 contain a nonzero number of points. 

4. Results 

4.1. Luminosity Evolution 

We apply the test of independence outlines in the § 3.1 to the entire BNB04 GRB sam- 
ple to test for luminosity evolution. For this analysis we use the flux threshold suggested by 
BNB04 of fmin = 0.5 photons cm -2 s -1 , decreasing the sample size to 985 bursts. Apply- 
ing this method, we find evidence for a strong 11.63 a correlation between luminosity and 
redshift. This evolution is well parameterized by a power law of the form g(z) = (1 + z) a , 
with an optimal value for the power law index (i.e when r(a)=0 given the transformation 
L — > V = L/g(z)) of a— 1.1 ± 0.3. The error estimates on a correspond to the 1 a range for 
this parameter given by the condition \r\ < 1. A plot of r(a) vs. a with the corresponding 1 
a levels are shown in Figure 6. These findings indicate that the average luminosity (modulo 
a beaming factor d£l) of GRBs in the universe has evolved with time. Because of the lack 
of beaming information, it may also be possible that the luminosity is remaining constant 
while the beaming factor dQ is actually evolving. As will be discussed in § 5, there is no 
observational evidence to suggest that this is the case. 

It should also be noted that r(a) appears to be strongly affected by the choice of the 
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flux threshold assumed for the sample. Plotted in Figure 6 is the optimal value for a vs. 
fmin- Not surprisingly, if we assume no flux threshold (i.e. f min = 0), r approaches the 
overestimated value received from the standard Kendell r statistic, a similarly approached 
the value obtained by simply performing a power-law fit to the truncated data, a decreases 
steeply with increasing f m in, never reaching a stable plateau as one would hope would happen 
as the fmin approaches the true threshold of the detector. This underscores the importance 
of having a good understanding the thresholds of the detector used to collect the sample. 
BNB04 make a strong case for a threshold of f m i n —0.5 photons cm~ 2 s _1 based on where 
they see a strong drop off of detected events in the L — Z plane (see their Figures 6 & 6) 
and we adopt this value for all analysis presented in this paper. 



4.2. Luminosity Function 

The deduced parametric form describing the luminosity evolution allows us to use the 
C-method on the uncorrelated parameters L' and Z to obtain the cumulative luminosity 
function <&(!/). Shown in Figure 6 is the cumulative <&(!/) distribution plotted as $(> V) 
as a function of L' for all 985 bursts. Because the luminosity evolution has been explicitly 
removed, this distribution represents the luminosity function in the present epoch. Fitting 
a double power law to the curve yields $(> V) oc L'~ 0S23 and $(> V) oc L'~ 19m for 
the low and high luminosity ranges respectively, separated by a break at a luminosity of 
roughly ~ 10 52 . These slopes are very similar to those reported by Lloyd-Ronning, Fryer, 
& Ramirez-Ruiz (2002) who found a GRB cumulative luminosity function with power law 
slopes of ki = —0.51 and k 2 = —2.33 below and above a break at about V = 5.9 x 10 51 . 
These values can also be compared to the luminosity functions found by Maloney & Petrosian 
(1999) who employ the C-method to account for selection effects in quasar samples. They 
find that the quasar luminosity function exhibits a double power law form with indices of 
k\ = —1.16 and k 2 = —3.59. 

Next, we differentiate the cumulative luminosity function with a 3-point Lagrangian in- 
terpolation to find the differential luminosity function d§/dL', or what is commonly referred 
to as simply the luminosity function ip(L'). This function represents the total number of 
bursts with luminosity between V and V + dV . A plot of the ip{L') vs. L' is shown in 
Figure 6. The function falls roughly as ip(L') oc L'~ L5 below the break energy of ~ 10 52 
with a sharp decline for higher L' . This power law index is identical to the slope found by 
Lloyd-Ronning, Fryer, & Ramirez-Ruiz (2002) who found L'~ L5 and similar to the index 
found by Schaefer (Deng) who found / y - 1 - 7±01 from (L, z) data estimated from a combined 
use of the lag-luminosity function and variability-luminosity function, although the latter 
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did not account for any selection biases in their data set. This value is also similar to results 
of several studies that used the measured flux distribution with an assumed density distri- 
bution p(z), such as (Schmidt 2001) who uses the star formation rate to infer a p(z) and 
finds ip(L') oc L~ 1A . The shape of the GRB luminosity function has important implications 
to jet model theories which predict specific power law indices for various jet structures. A 
comparison between theorized shapes and our deduced values will be discussed in more detail 
in §5. 



4.3. Density Evolution 

Using the alternative definition of the associated set, we can construct the cumulative 
density distribution a(z) = f Q z p(z)(dV/dz)dz, or the total number of GRBs per comoving 
volume, up to a given redshift. The cumulative distribution is shown in Figure 6 plotted as 
er(> z) as a function of z. The distribution of GRBs appears to increase smoothly with z, 
without a pronounced break at any distance, but with a flattening at high redshift indicating 
a drop off of events between 5 < z < 10. To get a better look at this density evolution, we 
can plot the cumulative density distribution a(z) as a function of comoving volume V(z) as 
seen in Figure 6 If the density of GRBs per comoving volume V(z) is constant, i.e. p(z) = p , 
then it should follow that a{z) oc V(z). We can test for evolution by fitting a{z) vs V(z) 
to a simple power law a(z) oc V(z)P and looking for deviations from the constant density 
case. An index of f3 ^ 1 indicates the presence of density evolution, with f3 > 1 and (3 < 1 
signifying an increasing and decreasing population respectively. Using the definition of V(z) 
in a flat universe of 
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h j ^/n m (i + z ) 3 + n A _ 

we find that the cumulative density distribution increase with z roughly as a(z) oc U 125 at 
low redshifts before falling off at higher redshifts. From these results we can deduce that the 
GRB density has undergone complicated evolution, increasing as p ~ \/ - 25 before peaking 
between z ~ 1 — 2 and then decreasing. To obtain a more quantitative look at the shape 
of the comoving rate density p(z), we again use a 3-point Lagrangian interpolation routine 
on cr(z) to find the differential cumulative distribution da/dZ. We can then convert this 
differential distribution into a comoving rate density through the relation: 

/ ~n da s / dV\ 1 „ s 

In Figure 6 we show the resulting comoving rate density plotted as a function of z. It can be 
seen that the GRBs density function increases out to a redshift between 1 < z < then flattens 
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before beginning to show signs of a turn over at a redshift of z > 3. This is in contrast to 
previous estimates of the comoving rate density by Schaefer (Deng), Lloyd- Ronning, Fryer, 
& Ramirez- Ruiz (2002), and (Yonetoku et al. 2004) all of who find a flattening of the GRB 
population with no apparent turn over out to a redshift of z ~ 10. It is also in contrast to 
results reported by Murakami et al. (2003) who also used the lag-luminosity correlation to 
estimate the GRB formation rate. There the authors find the GRB formation rate increases 
steadily out to a redshift of at least 4, but it should be noted that this work did not take 
into account the detector selection effects discussed above so a direct comparison may not 
be appropriate. As opposed to these previous findings, the turn over observed in our data 
quantitatively matches the global behavior of the star formation rate of the universe which 
has been observed to peak between 1 < z < 2 followed by a steady decline (Madau et al. 
1996; Steidel et al. 1999). A more detailed comparison between the GRB comoving rate 
density and the supernova and star formations rates will be continued in §5. 

5. Discussion 

We find an 11.63 a correlation between the luminosity and redshift data deduced from 
the lag- luminosity correlation, strongly suggesting an evolution of the average luminosities 
of GRBs. We show that this correlation can be parameterized as a power law as L(z) = 
(1 + z) L7±a3 . This value agrees extremely well with the results presented in Lloyd-Ronning, 
Fryer, & Ramirez- Ruiz (2002) who found a power law index of a — 1.4 after performing a 
similar analysis on (L, z) data estimated using the variability-luminosity correlation. These 
results imply that the average energy emitted per unit time per unit solid angle by GRBs 
was much higher in the distant past compared to relatively recent events. This is consistent 
with previous studies showing strong source evolution and also recent observations of under 
luminous nearby GRBs. Due to our lack of knowledge regarding the beaming angle of the 
bursts in our sample, it is also possible that the increase in the apparent luminosity is due 
to an increasing collimation at higher redshifts. As we will discuss in more detail below, we 
find no evidence for beaming angle evolution in the current sample of GRBs with known 
redshift and jet opening angle, suggesting that this increase in luminosity can not be due 
simply to an evolution of the collimation of the gamma-ray emission. 

5.1. Comparison to Other Objects 

Such a steep luminosity evolution is not uncommon in other astrophysical objects that 
show evolution with redshift. Maloney & Petrosian (1999) perform a similar analysis using 
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the statistical techniques described in this paper on a combination of several quasar samples 
and find that the quasar luminosity function evolves as L(z) = (1+z) 2 - 58 up to a redshift of at 
least 2. There is evidence that this evolution may then become constant up to a redshift of at 
least 3 (Boyle et al. 1993). We find no such break in the luminosity evolution of GRBs, which 
in our case can be adequately fit by a single power law between at least < z < 10. The 
authors also find a density evolution of a(z) oc V 1A9 similar to the power law of a(z) oc V 1 ' 25 
that we find in GRBs at low redshifts. A more detailed look at their comoving rate density 
estimate shows that the quasar density rises as p ~ (1 + z) 2 - 5 before peaking at z rs 2 and 
then declining rapidly as roughly p(z) (1 + z)~ 5 for z > 2.0. This is qualitatively similar to 
the trend we deduce from the GRB sample, which rises as p ~ (1 + z) 2A to a z 1 although 
the proceeding decline is much more shallow as p ~ (1 + z)~ os and extends to at least a 
redshift of z ~ 6 before dropping off sharply. 

There is also evidence for significant evolution in the luminosities of star forming galax- 
ies, which is perhaps a more relevant comparison to GRBs because of their suggested as- 
sociation with active star forming regions (Djorgovski et al. 1998). Hopkins (2004) used a 
compilation of recent star formation rate density measurements as a function of redshift to 
constrain the evolving luminosity function of star-forming galaxies. He finds that the pre- 
ferred evolution in a standard cosmology is given by L(z) = (1 + z) 2 - 70±0 - 60 out to a redshift 
as high a z ~ 6. At the same time he finds evidence for a very shallow density evolution given 
by p(z) ~ (1 + z )°- 15±om , markedly different from steep density evolution p(z) ~ (1 + z) 2 - 5 
that we estimate for GRBs at low redshift. This would indicate that GRB luminosities have 
evolved at a slower rate, but that their density in the past rises much more steeply compared 
to the number of star forming galaxies. It could also mean that the number of GRBs per 
star forming galaxy has evolved rapidly with cosmic time. 

Perhaps more interesting is the comparison between the GRB comoving rate and lumi- 
nosity densities with the global star formation rate history of the universe. Because GRBs 
suffer little extinction and are potentially detectable out to redshifts of z ~ 10, they could 
offer a unique tracer to the SFR history. They would allow for a more complete sampling 
of dust enshrouded star forming regions that may be missed in traditional SFR estimates 
based on the UV "drop-out" technique that is currently employed to identify Lyman break 
galaxies. The shape of the SFR at low redshifts z < 1 is relatively well understood, showing 
an order of magnitude increase from < z < 1 (Madau et al. 1996; Fall, Chariot, & Pei 
1996). These early estimates suggest that the star formation activity peaks around z ~ 1 
followed by a rapid decline at higher redshifts. However, further observations of hundreds 
of Lyman break galaxies at redshifts of z ~ 3 and 4 have shown that the SFR may remain 
constant after reaching a maximum around 1 < z < 2 (Steidel et al. 1999). Recent deep 
surveys with the Subaru (Iwata et al. 2003) and Hubble Space Telescopes (Bouwens et al. 
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2003) out to z ~ 5 and 6 show evidence for a mild evolution of the SFR at redshifts z > 3, 
with measurements based on photometric redshifts showing a constant SFR out to z ~ 6 
(Fontana et al. 2003). These recent SFR estimates qualitatively match the deduced GRB 
comoving rate density shown in Figure 6. At low redshifts z < 1, the GRB rate density 
increases as p(z) ~ (1 + z) 2 ' 5 roughly matching the rise in the SFR over the same range, 
with a peak somewhere between 1 < z < 2. The following flattening and decline between 
2 < z < 6 in the GRB p(z) matches the global properties of the SFR estimated from the 
recent deep surveys. 

Of course the comparisons between the GRB comoving rate density and the SFR are 
simply phenomenological, since we have as of yet no way of connecting the amount of star 
formation for a given amount of GRBs. Ultimately this conversion factor depends on knowl- 
edge of the GRB progenitor and the initial stellar mass function (IMF) and how it changes 
with redshift. In the case of the collapsar model (Woosley 2000; MacFadyen & Woosley 
1999), the rate of GRBs produced for a given SFR would increase sharply with redshift, 
as is the case for all core collapse events, due primarily to the redshift dependence on the 
IMF. However, the connection between the GRB p(z) and the SFR would be straightforward 
since the progenitors would consist of massive stars with short lifetimes making them direct 
indicators of the SFR at that redshift. If the mass range of the progenitors and the redshift 
dependence of the IMF is know (or assumed) then it would be possible to calculate a con- 
stant that directly relates p(z) to the SFR. The case is more complicated for binary merger 
models since there would be a substantial delay between the formation of the progenitor star 
and the final merger event that produces the GRB. The distribution in the delay times is 
not well known for SNe events, much less GRBs, but it is expected to be large enough to 
dissociate the GRB p(z) and the active SFR at a given redshift. The peak we observed in 
our deduced values for p(z) matching the peak of the current SFR estimates hardly seems 
coincidental, tentatively favoring the core collapse models. 

It is interesting then to compare our demographic results to that of various types of 
supernovae. There is overwhelming observational evidence and theoretical discussion sug- 
gesting a GRBs-SNe connection, including observations of supernovae bumps in afterglow 
lightcurves (Stanek et al. 2003; Hjorth et al. 2003) and a deduced collimation corrected en- 
ergy that is narrowly clustered around the typical SNe energy of 10 51 ergs (Frail et al. 2001; 
Bloom 2003). Although the intrinsic luminosity of type la SNe are a priori assumed to be 
constant with redshift (hence no luminosity evolution), we can still compare the formation 
rates of SNe Ia/b/c to GRBs, although the b/c events are obviously of more relevance to 
GRB models. Unfortunately, very little SNe data is available for the high z universe with 
only 7 SNe at z > 1.25 of the 42 SNe detected in the redshift range of 0.2 < z < 1.6 by 
the Advanced Camera for Surveys (ACS) on the Hubble Space Telescope (Riess et al. 2004). 
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Data on core collapse supernova accounts for only 17 events of this sample, going out to a 
maximum range of z ~ 0.7. Dahlen et al. (2004) use this data to estimate the core collapse 
SNe (CC SNe) rate between 0.3 < z < 0.7 and find a steep (about a factor of ~ 7) increase in 
the SNe p(z m 0.7) compared to the local rate presented by Cappallaro et al. (1999). Shown 
in Figure 6 are their data points for CC SNe plotted over the GRB comoving density, both 
rates normalized to 1 at z — 0.7. The two data points, although limited, do agree with the 
rise of the GRB p(z). A direct SNe-GRB comparison at higher redshifts will have to wait 
until the launch of the SNAP spacecraft which is predicted to find thousands of supernovae, 
including a significant number at high redshift. 



5.2. The Nature of the Luminosity Evolution 

The observed luminosity evolution that we observed in the (L, z) data leads to the 
conclusion that the GRB progenitor population has most likely evolved in such a way as to 
create more energetic or more narrowly beamed bursts in the distant past. Speculations on 
the nature of this evolution are dependant on the progenitor model and how the properties 
of their population are affected by the conditions of the early universe. In the case of highly 
rotating massive stars (i.e. the collapsar model), the overall progenitor mass and/or rotation 
rate could be the determining factor. There is ample evidence suggesting that the so called 
population III stars were much more massive than their present day counterparts. This is 
suggested by recent work showing that the stellar initial mass function (IMF) has evolved 
with time, having a much higher value in the distant past. This higher IMF is due to 
various factors, although it primarily is due to the lower metallicity in the early universe. 
The amount of material lost to stellar ejecta has also been shown to be dependant on the 
stellar metallicity, causing these early stars to retain more of their mass until their eventual 
collapse. 

Although the relationship between progenitor mass and emitted energy and/or beaming 
angle is not straight forward, there are reasons to think that this increase in average mass 
could result in an increase in the total energy budget available to a burst. MacFadyen & 
Woosley (1999) show that under the right conditions, the collapsar model could produce 
more energetic bursts with increasing stellar mass, up to some limit dictated by the energy 
needed by the GRB jet to punch through the stellar envelope. Proponents of black hole 
accretion disk models have also shown that the rate of accretion onto the central engine of 
the GRB increases dramatically as a function of the progenitor mass, increasing the overall 
energy available to the burst. 

Unfortunately, a simple increase in the overall energy budget cannot by itself explain 
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the deduced luminosity evolution. Frail et al. (2001) and Bloom et al. (2003) have recently 
shown observational evidence suggesting that the collimation corrected GRB energies E 1 
are actually narrowly cluster around the 10 51 ergs typical of SNe explosions. They come to 
this conclusion by correcting the observed prompt isotropic equivalent energy release E iso 
of several GRBs with known redshift by a factor of 1 — cos 9j, where 9j is the canonical jet 
opening angle. These angles are derived from broadband breaks observed in the afterglow 
light curves attributed to the slowing of the GRB jet to the point where the relativistic 
beaming angle of the radiation 1/T becomes greater then 9j. Bloom et al. (2003), using a 
larger sample of bursts, show that the corrected energies cluster around 1.33 x 10 51 ergs with 
a variance of 0.35 dex, or a factor of 2.2. Guetta et al. (2003) has reported a similar result 
when correcting for the isotropic luminosity Li so , although not as narrow as the E 1 results. If 
the collimation corrected energy and luminosity are indeed invariant with redshift, it directly 
implies that the brightening of the apparent isotropic equivalent luminosity is actually due 
largely to an increase in the beaming factor as a function of redshift and not an increase in 
the overall energy of the burst. There are physical arguments that can be made in the case 
of the collapsar model that would suggest that more massive progenitor stars could indeed 
produce more collimated jet outflows. 

Plotted in Figure 6 are the E 1 estimates from Friedman & Bloom (2005) for a little over 
two dozen GRBs. Furthermore, plotted in Figure 6 is the canonical jet opening angle for the 
same two dozen GRBs. By apply a standard Kendall rank order r statistic we can measure 
the degree of correlation in these two samples in a nonparametric way (i.e. without assuming 
an underlying correlation type). We find a correlation strength of r = 0.093 between E 1 and 
redshift and r = 0.163606 between 9j and redshift, where a tau of 1 signifies a significant 
correlation. Therefore, there is no deduced redshift dependency that would suggest any 
evolution of the jet opening angle or E 1 with redshift in the pre-Swift data set. This lack 
of redshift dependency stands to be tested in the Swift era as more GRBs with measured 
jet break times are observed over a broader redshift range, but if it is confirmed then it 
would imply that the evolution of some jet property other than the collimation factor must 
be responsible for the brightening of GRBs with redshift. Speculations on the nature of this 
evolution are dependant on the jet model used to explain the emission. The simplest model 
assumes a uniform energy distribution per solid angle e(9) across the jet with a sharp drop 
beyond 9j. In this scenario, the observed distribution in GRB energies is directly due to the 
diversity of jet opening angles, as is the observed values of the jet break time tj. The lack 
of any concrete evidence for an evolution of 9j with redshift combined with the observation 
that the collimation corrected E y and L 1 are very narrowly clustered, create difficulty for 
the uniform jet model to explain any kind of evolution in luminosity. One of the observed 
conditions above would have to be broken in order to accommodate any such evolution with 
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this model. 

In a structured jet model, the GRB jets are identical having a quasi-universal shape with 
a fixed opening angle and a nonuniform energy distribution per solid angle. The diversity 
in the observed jet break time and isotropic energies would then be a result of varying 
viewing angles away from the jet axis V . Furthermore, an observer viewing the GRB at a 
small V would see an extremely powerful burst, with the observed luminosity declining as 
some function of increasing V . A jet structure with a functional form of e{0)~ 2 is required 
to reproduce the observed tj oc E iso , i.e the Frail et al. (2001) and Bloom et al. (2003) 
results. If the requirement of a narrow E 1 and L 7 distribution is broken, then any power 
law structure 9 k could still produce the observed steepening in the afterglow light curve. 
In this case, the luminosity evolution would manifest itself not as an narrowing of 9j but 
rather as an overall increase in the normalization of e(0j). Another possibility would be an 
evolution of the morphology of e(theta) as a function of redshift. If the power law index 
e(6j) ~ 9j has evolved with time, or if e(0j) has evolved from a non-power-law shape (e.g. a 
Gaussian profile), such that e(0j) varies more slowly with viewing angle, then there would 
be a markedly different luminosity distribution at high redshift. A third, rather implausible, 
explanation is a preferentially small viewing angle V at higher redshift, although there is no 
physical reason to think that this is at all possible. Therefore, it would seem that evidence 
of luminosity evolution in the presence of the observation that E 1 and L 7 are narrowly 
distributed and the lack of any evidence of an evolution of 0j with redshift, favors a quasi- 
universal jet model over a uniform jet model. This is primarily due to the inability of the 
uniform jet model to explain any kind of luminosity evolution with redshift without a parallel 
evolution in the jet opening angle, something that is not currently observed. 



5.3. Luminosity Functions and Jet Model Discrimination 

Because the energy distribution r)(0j) of the structured jet model is well defined, it 
can make specific predictions regarding the GRB luminosity function 0(L). In the case of 
power-law structured jets t{9j) oc 9j k , resulting in a predicted luminosity function with a 
slope of7 = l + 2//c (Zhang & Meszaros 2002). The "canonical" k = 2 model would yield a 
luminosity function oc L~ 2 (Rossi, Lazzati, & Rees 2002), whereas the quasi-universal gaus- 
sian structured jet model predicts oc L^ 1 (Lloyd- Ronning, Dai, & Zhang 2004). Although 
the uniform jet model also exhibits a well defined e(9j), it cannot make any firm predictions 
about the shape of the GRB luminosity function due to the random variation 9j. 

The luminosity function deduced from our analysis is well fit by a single power low 
with an index of oc for luminosities below about 10 52 , with a sharp decline for higher 
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luminosities. This value is intermediate between the expected value for the k — 2 power law 
and gaussian structure models. An index of L'~ L5 actually predicts a power law model with 
k = 4 which is much steeper than the k = 2 shape needed to keep E iso 9 2 ~ constant. It is also 
possible that a simple gaussian or power law profile for rj(9) is simply an oversimplification. 
It has been pointed out that by Lamb & Graziani (2003) that r](9) would have to fall off 
steeper than 9~ 2 at large angles if the quasi-uniform jet model is to explain the relative 
numbers of x-ray flashes to GRBs. This may explain why so many studies have found a 
sharp decline above some break energy possibly indicating a different value for k at low and 
high luminosities, i.e large and small opening angles respectively. In any case, it would not 
be unfeasible to think that the jet morphology is more complicated than a simple power or 
gaussian profile, as is suggested by our results. 



6. Conclusions 

In this work we perform demographic studies on a large sample of luminosity and red- 
shift data found through the use of the lag-luminosity correlation. By applying maximum 
likelihood techniques, we are able to obtain an estimate of the luminosity evolution, lu- 
minosity function, and density distributions in a way that accounts for detector selection 
effects. We find that there exists a strong (11.63 a) correlation between luminosity and 
redshift that can be parameterized as L(z) = (1 + z) 2 ' 58 . The resulting cumulative $(-£/) 
and differential dQ/dU luminosity functions are well fit by double power laws separated by 
a break energy of about 10 52 ergs s _1 , with d§/dL' exhibiting a power law shape of L~ L5 
below this luminosity. This value does not immediately discriminate against any proposed 
structured jet models, but it may indicate that a more complicated jet profile is need to 
explain the observed luminosity function of GRBs. The GRB comoving rate density is found 
to increase as p-y(z) oc (1 + z) 2 ' 5 to a redshift of z ~ 1 followed by flattening and even- 
tual decline above z > 6. Although, the conversion between p 1 {z) and an estimate of the 
SFR cannot be quantitatively made due to the uncertainty regarding the GRB progenitors 
and their initial stellar mass functions, it can be said that p-,(z) does qualitatively follow 
recent photometric estimates of the SFR, as would be expected from massive short lived 
progenitors. We stress that these conclusions are based on the validity of the lag-luminosity 
correlation, which still stands to be confirmed as new redshift data becomes available. A 
full confirmation, and most probably further calibration, of this distance indicator will have 
to wait for addition detections with the Swift spacecraft, which should have the collecting 
area necessary to obtain high signal-to-noise energy dependant light curves from which to 
measure statistically significant lags. Even with the slew of directly measured luminosity 
and redshift data expected to come from the Swift mission, empirical distance indicators 
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may still play an important roll in expanding the available GRB data set through the use of 
archival BATSE and BepooSAX data for statistical studies such as the analysis performed 
in this work. 
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Figure Captions 



Fig. 1. - The luminosity and redshift data used in our analysis as deduced from the 
lag-luminosity correlation. The dashed line represents an imposed 0.5 photons cm~ 2 s _1 . . . 
cut to the original 1438 bursts analyzed by BNB04, producing the sharp cutoff in the data. 
This leaves a total of 985 bursts with a median redshift of 1.64. 

Fig. 2. - A plot of the lag-luminosity plane for the events under consideration. The error 
in the lag and flux measurements are used to determine the uncertainty in the luminosity 
values which in effect is used as a weight in the maximum likelihood technique that estimates 
the correlation strength between L and z. 

Fig. 3. - A representation of the associated sets used in the Lynden-Bell technique. For 
each data point with (Lj, zi), the solid line represents the minimum luminosity or maximum 
redshift that the burst could have had and still have been observed. Employing the Lynden- 
Bell technique with associated sets defined as Nj (Lj <L<oo,0<z< z rnax (Li)) produces 
the cumulative distribution for the vertical axis, whereas the Mj associated set produces the 
distribution for the data represented on the horizontal axis. 



This preprint was prepared with the AAS IATj^X macros v5.2. 
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Fig. 4. - Generated luminosity and redshift data used to test the ability of the Efron 
& Petrosian method to estimate the correlation strength and functional dependence for data 
of a given flux cut. 

Fig. 5. - A histogram of the difference between the known correlation index and the 
reconstructed index (p — q) for a large set of generated (L, z) data with arbitrarily imposed 
flux thresholds. The error in the method is tightly centered around p — q = 0. 

Fig. 6. - The correlation statistic r is plotted as a function of power law index 
g a {z) = (1 + z) a parameterizing the luminosity evolution. The solid line represents the 
a index that minimizes the correlation between L' and z. The dotted lines show the 1 a 
statistical error in the a parameter. 

Fig. 7. - The correlation parameter a is plotted as a function of the flux cut applied 
to the BNB04 sample. The optimal a is highly dependent on the choice of the cut value, 
showing the importance of a good understanding of the detector flux threshold. 

Fig. 8. - The luminosity and redshift data used in our analysis as deduced from the 
lag-luminosity correlation. The dashed line represents an imposed 0.5 photons cm~ 2 s _1 . . . 
cut to the original 1438 bursts analyzed by BNB04, producing the sharp cutoff in the data. 
This leaves a total of 985 bursts with a median redshift of 1.64. 

Fig. 9. - The present epoch GRB luminosity function tp(L') = d$/dL', representing 
the number of events between the luminosity V and V + dV . The function falls as roughly 
oc I/ -1 - 5 for luminosities below the break energy of 10 52 ergs s -1 . 

Fig. 10. - The cumulative density distribution a(z) = J* p(z)(dV/dz)dz, representing 
the total number of events up to a given redshift. The flattening at high redshift indicates 
a drop off of events around z ~ 5 — 10. 

Fig. 11. - The cumulative density function a(z) plotted as a function of comoving 
volume. The dashed line represents the increase in the number of sources if the GRB density 
were constant throughout the history of the universe. The GRB density has increased as 
p ~ yo.25 k e f ore peaking between z ~ 1 — 2 and then decreasing. 

Fig. 12. - The comoving rate density p(z) as a function of redshift. The rate density of 
sources can be seen to follow the evolution deduced from Figure 6, increasing to a redshift of 
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1-2 then flattening before decreasing at higher redshifts. The circles represent high redshift 
cc rates from Dahlen et al. (2004) whereas the square point is the local rate found Cappallaro 
et al. (1999). The increase in the cc event rate with redshift qualitatively matches the overall 
increase in the GRB comoving rate density. 

Fig. 13. - Plot of the collimation corrected total emitted energy of 23 GRBs with 
known redshift and beaming angles. No significant correlation can be seen with redshift. 

Fig. 14. - The beaming angles 9j of 23 GRBs with known redshift. The lack of a 
significant correlation with redshift is quite evident. 
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